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Abstract 

This paper is devoted to the numerical resolution of an anisotropic 
non-linear diffusion problem involving a small parameter e, defined as the 
anisotropy strength reciprocal. In this work, the anisotropy is carried 
by a variable vector function b. The equation being supplemented with 
Neumann boundary conditions, the limit e— >0 is demonstrated to be a 
singular perturbation of the original diffusion equation. To address effi- 
ciently this problem, an Asymptotic-Preserving scheme is derived. This 
numerical method does not require the use of coordinates adapted to the 
anisotropy direction and exhibits an accuracy as well as a computational 
cost independent of the anisotropy strength. 
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1 Introduction 



The aim of this paper is to build an efficient numerical method for solving an 
anisotropic diffusion problem where the anisotropy is carried by a vector b. 
This work is motivated by investigations of strongly magnetized plasmas, more 
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specifically the study of the Euler-Lorentz model in a low Mach number regime 
and in the presence of a large magnetic field. This framework is characteristic 
of the magnetically confined plasma fusion [T7J [551 In this context, the 

asymptotic parameter e represents the gyro-period of particles as well as the 
square root of the Mach number, the vector field b being the magnetic field 
direction. Therefore the e values can be very small in some sub-regions of the 
computational domain where the magnetic field is large, inducing then a severe 
anisotropy of the medium, while being large in other sub-domains for interme- 
diate and small strength of the magnetic field. Another important property of 
this system is the time dependence of the magnetic field defining the anisotropy 
direction. These two main characteristics define the framework of the present 
paper whose purpose is to design a numerical scheme for anisotropy ratios rang- 
ing from eCl to e~ 0(1) and for a time varying anisotropy direction. In order 
to address efficiently these requirements, the numerical method should not rely 
on a coordinate system adapted to this anisotropy direction. The use of adapted 
coordinates would imply mesh modifications accordingly to the evolution of b, 
an intricate and expensive procedure we wish to avoid. Thus, the numerical 
method introduced here will carry out the anisotropic non-linear diffusion prob- 
lem on a mesh independent of the anisotropy direction. 

This scheme will be detailed on the following model problem 



In this system, f2 is a bounded subset of M. d (d= 1,2,3), e > is a fixed constant 
parameter and, for any xG<9f2, u = u(x) stands for the unit outward normal 
vector. V x and V x - stand for the gradient and the divergence operators with 
respect to the space variable x. We assume that H e : f2 — > K*j_ , b:f2— >-M d — 
{0}, / e :Q— >R, S e :f2— s-R d are given and the unknown of the problem is the 
function p e : f2 — > R. The tensor product of two vectors u and v is denoted 
u®v. Finally, we assume that, for any e, the function pt-^g e (p) is strictly 
increasing and can be non-linear. This equation is well suited for the plasma 
fusion context above depicted. It allows the computation of the plasma pressure 
in order to guarantee that the forces vanish in the low Mach regime for strongly 
magnetized plasma. The function denoted g e defines the internal energy of 
the fluid with respect to the pressure. This relation may be non-linear, which 
motivates the investigation of non-linear anisotropic problems. However, the 
derivation of this equation is out of the scope of the present paper and we refer to 
related works (see [51 USE]) for detailed explanations. Furthermore, we wish to 
present the numerical method in a context wider than the strict plasma context, 
since anisotropic diffusion problem are encountered in many applications. Good 
examples of these applications are, for instance, image noise filtering, convection 
dominated diffusion equations and more generally diffusion problem with strong 
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medium anisotropies. The model equation (|1.1[> is representative of a large 
enough variety of problems, up to slight changes, and will be considered to 
detail the numerical method. 

Developing an efficient numerical method to compute the solution of this 
diffusion problem, regardless to e values, is a difficult task. Indeed, the limit 
e — > is a singular limit for the problem (|I.I[) , the diffusion equation degenerating 
into the following one 



The system (|1.2p is ill-posed, its solution being non-unique. More precisely, if p 
is a solution of (|I.2I) and c:f2— >-R is a function verifying b-V x c = on Q, then 
po + c defines a new solution of (II .2[) . However, p the limit of p e solution of 
(jl.ip is uniquely defined by the limit problem as demonstrated in Section [2J but 
a direct discretization of the diffusion problem (|1.1[) gives rise to a linear system 
with a conditioning number that blows up for vanishing e. This property has 
been outlined for numerical studies of elliptic equation singular perturbations 



To tackle this difficulty, an Asymptotic-Preserving (AP) scheme is introduced 
to compute the solution of the anisotropic diffusion problem for e — 0(l) and to 
capture po, the solution of the limit problem, for small e values. This property 
should be provided without any limitations on the discretization parameters 
related to the value of e. These requirements are compliant with the properties 
of AP-schemes originally introduced in [30] and developed in [32] for diffusive 
regimes of transport equations. These techniques have received numerous ex- 
tensions to other singular perturbation problems: relaxation limits of kinetic 
plasma descriptions [21 [37] , quasi- neutral limit of fluid and kinetic plasma 
models [H [101 EU [HI [13l HH [1H [23], hydrodynamic low Mach number limit 
[24l [31] , radiative hydrodynamics [7j [8] , fluid and particle flows [9] and strongly 
magnetized plasmas as well as heterogeneous media [H O [HI [HI H7J HOI 121] • 

The Asymptotic-Preserving property of the presented method is obtained 
thanks to a decomposition of the solution, introduced in [20] and also used in 
[U[S][T5]- It consists of the following identity p e — TT e + q e , ir e being the solu- 
tion mean part, with respect to the anisotropy (b) direction, q t the fluctuating 
part. These two components verify Tr e £K and q c £K i - 1 K defining the func- 
tions constant along the b-direction, K 1 - the functions of zero mean value along 
b. This decomposition was first developed for meshes adapted to the anisotropy 
direction [U[5D], f° r which, the discretization of K is straightforward. A direct 
discretization of the sub-space K is, on the other side, much more intricate. 
This difficulty is overcome thanks to the introduction of a Lagrangian multiplier, 
in order to penalize the zero mean value property of the functions belonging to 
K . The method is extended in [T5] for computations with meshes indepen- 
dent of the anisotropy direction. This is achieved by introducing two more 
Lagrangian multipliers to discretize the sub-spaces. The size of the linear sys- 
tem providing the problem solution is then significantly enlarged. However, this 
drawback may be corrected thanks to a slightly different decomposition. In [16] 



{ 



-V x -( J ff o (b«)b)(V x po-So)) = 0, in Q, 
(i/ o (b(g)b)(V x po-So))-^ = 0, onSfi. 



(1.2) 



(see [2 [20]). 
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the solution is decomposed in two non-orthogonal parts which allows the defi- 
nition of two sub-spaces whose direct discretization is readily obtained without 
any Lagrangian multipliers. The size of the linear system obtained with this 
approach is considerably lowered compared to the previous method [15] , This 
method has been extended in [341 133] to non-linear diffusion equations. The 
path followed in the present paper still relies on the decomposition in K and 
K , However, the discretization of these sub-spaces is achieved using a differ- 
ential characterization, similar to the one introduced in [5j. This finally allows 
the computation of the solution thanks to a second-order problem for 7r e and a 
fourth-order problem for q f . This latter problem, in the framework of Neumann 
boundary conditions considered in this paper, can be recast into two elliptic 
problems. 

The method proposed finally reduces to the computation of three standard 
elliptic problems for which very efficient solvers can be used (for instance multi- 
grid solvers) . For the former approaches [HI [16] , the equations providing both 
components are not classical elliptic equations and the resolution of the linear 
system requires more sophisticated solvers. This complexity is resource de- 
manding and may be challenging for realistic three-dimensional computations. 
Finally this paper also presents an extension to non-linear reaction diffusion 
problems, a class of problems that has never been investigated in the previous 
works. 

The paper is organized as follows: in Section 2, the decomposition method- 
ology is presented. The linear case, i.e. with g e (p)(x) = G e (x)p(x) where 
G e : Q — > is a given function sequence, is first investigated: more precisely, we 
describe the decomposition procedure in the specific case where G £ is a strictly 
positive constant denoted A e , then we generalize this procedure to any func- 
tion G e :fl — >M*\_ by using well-chosen Sobolev spaces. Finally the non-linear 
problems are addressed by invoking Gummel's iterative algorithm. Section 3 
is devoted to presentation of the discretization. Finally, the efficiency of the 
numerical method is demonstrated in Section 4. 

2 Scale separation and solution decomposition 

In this section a scale separation is introduced to ensure the Asymptotic- 
Preserving property of the scheme. This is achieved by transforming the sin- 
gular perturbation problem (11.11) into an equivalent system for which the limit 
e — > is regular. For simplicity reasons, the linear case with constant G e is 
first considered for detailing the decomposition method. In this framework, the 
singular nature of the limit e — > is outlined and the limit problem, providing 
Po = lim e _i.oPej is stated. A development to linear cases with variable positive 
functions G e is then presented and finally, thanks to Gummel's iterative method 
|28j . the non-linear case is addressed by using a sequence of linear problems. 
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2.1 AP-scheme derivation for linear problems 
2.1.1 A simplified framework: constant G e 

We assume here that the given sequence (<? e ) e >o is of the form 

0e(p)(x) = A e p(x), 

where A e > is a known constant for any e > 0. Then the diffusion problem 
writes 



The limit solution pg of the singular perturbation problem (|2.1|1 verifies the limit 
problem 



This algebraic equation admits a unique solution under the assumption 



a requirement that must be fulfilled by the numerical method. To ensure this 
property, the methodology consists in using a decomposition similar to that 
of [3 O HH [HI [20]. The solution is decomposed into 7r c its mean part with 
respect to the anisotropy direction and q e the fluctuating part, which exhibits 
the property to have a zero mean value along the anisotropy direction. These 
two functions verify Tr e £K and q e £ K , K being the kernel of the elliptic 
operator defined by equation ll.2l These properties are capitalized on, to isolate 
in the problem 12. II the macro scale (providing 7r e ) from the micro scale (giving 
q e ) and thereby, build the Asymptotic-Preserving scheme. The main difficulty 
of the procedure lies in the characterization of the sub-spaces associated to the 
different scales. In [4j [13j [15l [20] the property of the functions populating K 
or K are imposed by a penalization technique. The methodology developed 
in this paper operates a similar decomposition on to K and K- 1 , but with a 
different characterization of these sub-spaces. Here, we shape the technique 
introduced in [5 j for a very specific framework, in order to discriminate the 
functions in K and thanks to differential properties, providing thus, an 
easy discretization. 

With this aim, we introduce the following Sobolev spaces: 




(2.1) 




(2.2) 



ff e (b(g)b)(V x Pe-S e ) = 0(e), 



(2.3) 



V={p£L 2 (fl) :b-V x p€L 2 (Q)} , 
W={q£L 2 (n):V x -(hq)£L 2 (n)} , 
W = {q£W: (bq)-u = 0on 90}, 



and we define K C V as 



K = {n£V :b-V x 7T = 0on n} 
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The goal is to reproduce the function decomposition into its mean and fluctuat- 
ing parts. The functions of K correspond to the mean part and the complemen- 
tary part is demonstrated to belong to K 1 - . This is the purpose of the following 
theorem: 

Theorem 2.1 We denote by V x -(bWo) the subspace of functions 6>eL 2 (f2) 
such that 

3X61^0, = V x -(bx), (2.4) 
and we equip it with the usual norm on L 2 (il). Then: 

• Wq equipped with the norm ||p||iv = ||V X • (bp)|| L2 is a Hilbert space, 

• Vx-(bWo) is a closed subspace in L 2 (fl), 

• K is a closed subspace in L 2 (f2). 

• We have the orthogonal decomposition 

L 2 (n) = K®K ± 1 with K ± = V x -(bW ). (2.5) 

The demonstration of this theorem will be omitted. It can be readily adapted 
from that of Theorem 2.1 from 5 r As a consequence of this theorem, the 
decomposition 

p e =Tr e + q e ,ir e eK,q e eK ± , (2.6) 

exists and is unique for any e>0. Therefore finding the particular solution 
Po which is exactly the limit of (p c ) c> o is equivalent to find ttq and go as the 
respective limits of (7r e ) e >o and (q e ) e> Q. Then, our goal is now to find some 
equations for 7r e and q e which are well-posed for any value of e, including e = 0. 
For this purpose, the decomposition (|2.6[) is introduced into (|2.1[) . yielding 



-V x • {H e (b® b) (V x g e - S e )) + e A e (tt 6 +q,) - e/ e , in O, 
(F £ (b®b)(V x g e -S e ))-iv = 0, on dVt. 

The variational formulation on V writes 

/ # e (b-V x g £ )(b-V x 60dx + eA e / (Tr e + q e )9dx 
Jn Jn 

= e I f e 9dx+ I H e (b-S t )(b-V x 9)dx, 



(2.7) 



(2.8) 



for any test function 6gV. 

In order to exhibit the equation providing ir e € K, the variational formulation 
is tested against 9 6 K giving 



/ (\ e n e -f e )9d X = 0, 
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which means that A e n e — f e 6 K 1 - for any e > 0. According to Theorem 12 .11 there 
exists a function h e € Wo such that 

A e7 r e -/ e = V x -(bft e ). (2.9) 

This equation furnishes a means of computation for n e . Firstly, applying the 
differential operator b • V x onto (|2.9|) leads to an equation for h e 

-b-V x (V x -(b/i £ )) = b-V x / £ , in Q, 
{bh e )-u = 0, on 90, 1 ' 

then, 7r e is retrieved thanks to 

7T £ = ^[/e + V x -(bft e )]. (2.11) 

Note that the system (|2.10j) - (|2.11j) is well-posed and does not degenerate for 
any value of e > 0, including e = 0. It provides a means of computing the macro 
component of the solution regardless to e values. 

To derive an equation for q £ 6 K , we now assume that the test function 9 
in (|2.8p is in K 1 - . According to Theorem 12. 1[ there exist two functions x an d h 
in Wo such that 

9e = V x -(bg, (2.12) 

and 

= V x -(b X ). 

As a consequence, the variational formulation of (|2.7j) can be rewritten as fol- 
lows: 

/ H e (b- V x (V x • (bZ e ))) (b- V x (V x • (b X ))) dx + eA e / (V x • (bl e )) (V x • (b X )) dx 

= e // e (V x -(b X ))dx+ / J ff e (b-S e )(b-V x (V x -(bx)))dx. 
7n jo 



We recognize the variational formulation of 

A e b-V x (V x -(b/ e )) 

(2.13) 



f b • V x (V x • (H e (b® b) V x (V x • (bi e )))) - e A e b • V x (V x • (bZ e )) 

= -b-V x (e/ e -V x -(ff £ (b®b)S e )), inO, 
(if £ (b(8(b)Vx(Vx-(b/ e )))-i/ = (ir e (b®b)S 6 )-i/, on dfl, 

(bl e )-u = 0, ondn. 



Therefore, coupling this system with (|2.12|) . we recognize a complete definition 
of q e which is well-posed for any e>0, including e = 0. Moreover this compu- 
tation of q e is totally compliant with the condition (|2.3|) and guarantees the 
Asymptotic-Preserving property of the scheme. 

At this point, we have established a system of equations for 7r c and q t which is 
well-posed for any e > but also for e = 0. Then, solving the well-posed equations 
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(|2.10l) . (I2.13[) , (|2.11[) and (|2.12l) provides n and qo as the respective limits of 
(7r e ) e> o and (g e ) e >o when e— >0. As a consequence, the sum 7To + (?o is exactly 
the solution p of (12.21) . Furthermore, we can remark that the limit e— >0 is 
regular for the reformulated model (pl^ - (pl^ - (f23T|) - (f2~T2l) . 

2.1.2 Case with variable G £ 

In this paragraph, we extend the method we have presented to the general linear 
case, i.e. to cases where g t is of the form 

flr £ (p)(x) = G e (x)p(x). 

G £ :f2 — »K is given for any e>0, and is supposed to be strictly positive on f2. 
In such a case, the diffusion problem (|1.1[) writes 

-V x -(i? £ (b«)b)(V x p e -Se)) + eG £ p e = e/ e , in fi, 
( J H' e (b(8)b)(V x p e -S £ ))-^ = 0, onffi. 1 > 

The study of these cases is motivated by the fact that the use of Gummel's 
algorithm on the non-linear case leads to the resolution of a sequence of 
linearized problems which are similar to (|2.14j) . We refer to Section T2.2I for 
more details about the linearization procedure. 

In order to solve the linear problem (|2.14p for any value of e, we use the 
method presented in the previous paragraph. Firstly, we define L 2 (f2;G £ ) by 

L 2 (il; G e ) = jp : fi -> R, \\p\\ 2 L 2 {n . Gc) = jf G £ (x) |p(x) | 2 dx < +oc J . 

Then we introduce the following weighted Sobolev spaces: 

V t = {p e L 2 (fi; G £ ) : b • V x p e L 2 (O; G £ ) } , 
PF £ = {g £ i 2 (0;G £ ) : V x • (G £ b<z) e L 2 (r!;G £ ) } , 
Wo, £ = {<jeVF £ : (G £ bg)-^ = 0on 917}, 

and the set representing the functions constant along the magnetic field lines 

K e = {ire V t : b-V x 7r = 0on ft} . 

Following the methodology presented in the previous paragraph and in [5], we 
deduce 



L 2 (f2;G e ) 



IS 



Corollary 2.2 Wo, e equipped with the norm ||p||w ,« — ||V X • (G e bp)|| 
a Hilbert space and V x ■ (G £ b Wo,e) is a closed space in L 2 (i7;G £ ) . Furthermore, 
K e is also a closed space in L 2 (il;G £ ) and we have the orthogonal decomposition 

L 2 {tt-G t )=K t ®K^ with J^ = -^V x -(G £ bW^ £ ). (2.15) 

G £ 
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From the orthogonal decomposition (|2.15p . the solution p e of (|2.14p can be 
uniquely decomposed as 



p £ = 7r £ + q e , 7r £ € K e , q e <E K e 



(2.16) 



Then, if we identify the limits ttq and qo of the sequences (7T £ ) £ >o and (q £ ) £ >o, 
we will find the limit po of (p £ ) £ >o by taking po = 7r o + ( Zo- 

In order to identify a set of equations satisfied by 7r £ and q e , we follow the 
same procedure as in the previous paragraph: we multiply (12. 14|) by a test 
function 9 £ V e and we integrate over ft. By choosing 9 in K e or in , we 
prove that 7r £ and q e are respectively of the form 

vr £ = [/ £ + Vx • (G £ b/i £ )] , q £ = iv x • (G e bZ e ) , (2.17) 

Gr £ Cr f 



where /i £ and Z £ are solutions of 

-b-Vxf ^-V x -(G e b/i e )j =b-V x | 



(G £ b/i £ )-i/ = 0, 



and 



in f2, 
on dfl, 



(2.18) 



b-V, 



G, 



■Vx- (ff £ (b®b)V x ( 7^V x -(G £ bZ £ ) 



G f 



-ebV x l — V x (G e b/ e 



-b-V x — (e/ e -V x -(ff 6 (b®b)S e )) 



ff £ (b®b)V x I 

[ (G £ b; £ )-^ = o, 



G, 



■v x -(G £ b^ £ ) 



in f2, 



•i/=(ff e (b<g>b)S 6 )-i/,on 50, 
on 



(2.19) 



As in the previous paragraph, we observe that the equations (|2.17p - (|2.18p - (|2.19p 
remain well-posed for any e > 0. As a consequence, the particular solution po of 
the limit problem we are looking for is exactly the sum ttq + qo where ttq and qo 
are computed by solving (|2~T7| - ([2~T5| - (|2~TO| with e = 0. 

Furthermore, the resolution of the fourth order problem (|2 . 19[) can be re- 
placed by the successive resolution of two homogeneous Dirichlet type problems 
which are 



-b-Vx^V x -(ff £ bL £ )j+6L £ 
(H e bL e )-v = 0, 



eb- 



in CI, 



on Oft, 



(2.20) 
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and 

-b-V x ^V x -(G f b/ r )j =Z £ -b-S £ ,infi, (921) 
(G e bZ e )-i/ = 0, ondfl. 

2.2 AP-scheme derivation for non-linear problems 

Finally, we consider the general model fll.l[) given in the introduction when the 
function p^g t (p) is non-linear. When e goes to 0, the model becomes 

/ -V x • (H (b®b) (V x p - S )) = 0, in Q, 
\(i/o(b®b)(V x p o -So))-^ = 0, on 90. [ ' 

Due to the non-linearity of the function p t-> g e (p) the orthogonal decomposition 
method cannot be used. Then we choose to linearize the diffusion equation (11.11) 
by using Gummel's algorithm developed in [28]. This iterative method consists 
in the approximation of the solution p e by a sequence (p e ,N)N>o defined by 

Pe,N+l =Pt,N + 5e,N , (2.23) 

and initialized with an arbitrary p e ,o- In this method, each £) £i at is viewed as a 
small correction of p e ,N in order to obtain p e ^ + \. Then, assuming that p e ,iv-t-i 
is a solution of (jl.l[) . it holds that 

-V x - Lff e (b®b) ; 1 



e / 

>;2 



+.9 £ (p e ,w) + ^,Ar3 £ (p e ,Jv) + 0(<5 £ 2 JV ) = / e , in SI, (2.24) 

rr . . V x p £! Ar + V x 5 £ .Ar-S £ \ 

7J £ (b®b) ! ■ ^ = 0, on dfl. 



Then, neglecting second order terms in S e jy, we obtain a linear diffusion problem 
for <5 £ jv which writes 



(2.25) 



- V x • (H e (b <g b) (V x< 5 £ ,Ar - S £ .at)) + eG £ .Ar £ £iA r = ef €;N) in Q, 
(ff £ (b®b)(V x £ £ ,AT-S £!A r))-iv = 0, on aO, 

where G £j at, f e .N and S £i at are defined by 

Ge,N = g'<:(Pe,N), /e.JV = /e -ffe(?>e,iv) , S e ,iV = S £ - V x p £i jv • 

For each value of AT, the problem (|2.25|) is of the same kind as ()2.14|) . So we 
can solve it by applying the method described in the paragraph 2.1.2. 

This sequence of linearized problems can also be obtained from Newton's 
iterative method to solve 

F e (p e ) = 0, (2.26) 
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where the differential operator F e is defined as 

F e (p) = - V x • \H e (b ® b) VxP ~ S ^ + ,g e (p) - / 6 . 

Indeed, Newton's method for solving (|2.26l) writes 

DF t {p^ N )(p^ N+1 -Pc,n) = -F e (p e ,N) , 

where DF e (p) is the derivative in p of the differential operator F e (p) and is of 
the form 

DF e (p) (5) = - V x • ^i? £ (b ® b) + g' e (p) x J . 

3 Numerical method 

In this section, we present a numerical method which allows to solve the diffu- 
sion problems (j2.14j) and (ll.l[) by using the decomposition approaches we have 
presented. First, we introduce some notations which will be used for the con- 
struction of the scheme, then we present the scheme itself for the general linear 
case (|2.14[) . Finally, we present the discretized version of Gummel's algorithm 
for the non-linear case. 

3.1 Notations and definitions 

We consider a uniform mesh (xi,yj) defined by 

. . ■ a a Xjxiax Xmin A Umax Vmin 

Xi = x min +iAx, y j =y min +jAy,Ax= — — — , Ay= — — — — , 

and we assume that the simulation domain is = [x_i/2 ) a;jv a ,+i/2] x 
[y-i/2, VNy+l/z]- We also consider the following subsets of Z 2 : 

I={0,...,N x }x{0,...,N y }, 

7={-i,...,jv x +i}x{-i,...,jv w + i}, 

/, = {0,...,JV B -l}x{0,...,JV l ,-l}, 
Z={-l,...,N x }x{-l,...,N y }, 

and we consider the notation 

h = m.ax(Ax,Ay). (3.1) 

Since the decomposition method we have presented in paragraph 2.2 is based 
on variational formulations of the diffusion problem for p e and uses the duality 
between the operators pi— >h-V x p and pi— > V x - (bp), we choose to approach 
these differential operators by dh and dh,* respectively such that the duality 
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property is preserved at the discrete level. For this purpose, we define dh and 
dh * such that 



(dh,6)i+l/2,j+l/2 — b;+l/2,j + l/2 ' 



' #i+l,3+l + \ 



2Ax 

0«+l,j+l — + ' 



2Ay 



(3.2) 



for all = {6i. 



3 



and 



{dh,*X)i,j= X! 

a£{±l} 



(bxX)i+l/2,j+a/2 — (bxX)i-l/2,j+a/2 



2Aa 



+ £ 

a£{il} 



(&SlX)i+a/2,3+l/2 — (&j/X)i+a/2,j-l/2 

2Ay 



(3.3) 



for all X=(Xi+i/2,i+i/2)(ij )e T;- 
3.2 Linear problems 

We assume that the function pt-^-g e (p) is given by 

9e(p)(x,y) =G e (x,y)p(x,y) , 

where G e is analytically known. We also assume that the functions b, H e , f e 
and S e are analytically known and we consider the following notations: 

fe,i,j — fe{Xi,yj) , 
G e^i.j G e (Xj , yj ) , 
^6,4+1/2,3+1/2 = G t (x i+l /2,yj+l/2) i 
#6,1+1/2,3+1/2 = H e (%i+l/2, 2/i+l/2) ! 

S e ,i+i/2j"+i/2 = S e (a;i-|_i/2, 2/3+1/2) ! 
bj+i/2,3+1/2 = b(a;,+i/2, 2/3+1/2)- 



Then, the diffusion problem (|2.14l) can be discretized under the following form 
j (-dh,*(H e (d h p e! a P p-h-S e ) + eG e p e>a pp) i:j = ef et ij, V(i,j)el, 



{H e (d h p £ , app - b • S e ) (b ■ f)) i+1/2 J+1/2 = , 



(3.4) 



and the approximation of p e , a pp of p e is computed at the points (xi,yj) £ f2. 

Since dh and dh,* have been chosen to be dual operators, we follow the 
decomposition approach we have presented in Section 2.1.2 at a discrete level 
by using some discrete variational formulations of (I3.4p . Writing 



Pe,app,i,j — 7fe,app,i,j "T Qe,app,i,j , 
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with 7T ej0 p P satisfying 

(dh^e,app)i+l/2,j+l/2=0, V(i,j) G 7* 

n e ,app and q e ,ap P are completely denned by 

1 



'e,app,i,j - 



fe,i,j + (dh,*(G t /l e ,app))j j 



(3.5) 



Je,app,i,y ■ 



Ge,i,j 



{dh,*(G ele.app)) ; t j i 



where h e ^ a pp — {h t app ^1/2^+1/2)^ ^ e j^ and le,app — (le,app,i+i/2,j+i/2)^ij^ e j^ 
are computed by inverting the following systems: 



- I dh I -^-dh.*{G t h,:.app) 



i+l/2,j+l/2 
A 

G, 



U. 



e,app,i+l/2,j+l/2 — 0, 



and 



= -e 

- i£,app,i+l/2,j + l/2 =0, 
/ / 1 

- I ^/l I -^-dh,*(G e l ej ap P ) 



b-S e 



J+1/2J+1/2 



i+l/2,j+l/2 



j+l/2,j+l/2 



V(i,j)G/„ 



(3.6) 



V(»,j)e/., 
V(i,j)eZ\h 



(3.7) 



I ^e,app, i+l/2j' + l/2 



/ / i+l/2,j+l/2 
= (L e ,app -b-S e ) i+1 / 2 ,j+l/2! V(i,j) €/* 



(3.8) 



= 0, 



V(i,j)eh\h. 



3.3 Non-linear problems 



In this paragraph, we detail the discretized version of Gummel's algorithm pre- 
sented in Section 2.2. In order to initialize the loop, we compute the following 
initial datas: 

Pe,o,i,j = P<ifi{xi,yj), V(i,j)eJ 1 
S e ,i+i/2,j+i/2 = S e (a; i+ i/2, 2/7+1/2), V(i,j)eh, 

fe,i,j = fe{Xi,Vj), 'i{i,j)el ± 

H e ,i+i/2,j+i/2 = H e {x l+1/ 2,y j+ i/2), V(i,j)eJ*. 



Then, the iV-th iteration of Gummel's algorithm is set as follows: 
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• Step 1: assuming that 

Pe,N,ij, V(i,j)e7, 

are known, we have 

(b-S e ,jv)i+i/2,j+i/2 = (b-S e — 9hp e) jv)»+i/2,j+i/2, V(i,j) Gl*, 

/e,JV,»J = fe,i,j-9e(Pe,N,i,j), V(i,j) G J, 

G £ ,JV,i,j = 9e(Pe,N,i,j), V(i,j)eJ, 
^£,^,1+1/2^+1/2 = fl , e (Pe,AT,i+l/2,j+l/2) , V(i, j) G J* , 

1 

where P € ,JV,i+l/2,j+l/2 = j(Pe,Ar,*+l,j+l +Pe,JV,*+l,j + f>e,iV,i,.j + l +Pe,JV,i,j)- 

• Step 2: we compute /i £ ,jv,»+i/2,j+i/2 an d ^,jv,i+i/2,j+i/2 for all 
by solving 

9ft [ — d h ,«,(G e ,N h e . N )) ) 

V V Ge ^ /A+1/2J+1/2 

V V G WA+l/2, J+ l/2 
> ^e,JV,i+l/2,j+l/2 = 0, V(i,j) G J*\/*, 

| — — d h ,*{H e L e , N ) j +eL £ ,jv J 

V e ' W / / i+l/2,j+l/2 

= -e(dft(^J-b-S e J ,V(i,i)G/*, 

V V £ ' W / / J+1/2J+1/2 _ 

. ^€,JV,i+l/2,j+l/2 = 0, V(i,j) Gl*\I», 

and 

- ( 9ft | -^—d h ^ (G £;JV Z £;JV ) J J 

V V £,7V / / i+l/2,j+l/2 

= (-C'e,iV-b-S £ ,jv)i+l/2 ,j+l/2> Gi*, 

. ^,iv,i+i/2,j+i/2 = 0, V(i,j)€h\I*. 

• Step 3: we compute 5 e ^s,j for all G/ by using 

1 



Je,N,i,j — ' 



fe,N,i,j + (9ft,* (G £ ,AT (/l £ ,AT + h,N))) 



Ge,N,i,j 

and we obtain p €l jv+i,»j for all e/. 

• Step 4: we compute p £ ,jv+i,i,j for all by using the boundary 

condition 

(9ftp £ ,jv+i)i+i/2,j+i/2 - (b-S £ )j +1 / 2 ,j+i/2 = 0, V(i,j) G J*\I*. 
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4 Numerical investigations of the AP-scheme 

This section is devoted to numerical investigations of the Asymptotic-Preserving 
scheme derived in Sections [2] and [3] The validation procedure consists in manu- 
facturing p e , an analytic solution of the model problem which is compared 
to the numerical approximation p e , a pp carried out thanks to the AP-scheme. 
These experiments are performed in two dimensions using a uniform Cartesian 
mesh independent of the anisotropy direction. For simplicity purpose, the first 
numerical experiments are performed in the framework on the linear model, but 
the conclusions drawn from these investigations apply to the general non-linear 
problem. 

4.1 Numerical convergence of the scheme 

The first numerical tests aim at demonstrating the convergence of the AP- 
scheme regardless to the asymptotic parameter values. With this aim, an ana- 
lytic solution is manufactured for the problem in the linear case, i. e. with 
g € {j>) = G e (x)p(x). First, the expression for the anisotropy direction b and the 
functions G e and H e are defined on = [1,2] x [1,2] thanks to 

G e (x,y) = H e (x,y) = l + sin 2 (x) sin 2 (y) , (4.1) 
b= (sin#, — cos#) with &{x,y) = arctan I — I . (4.2) 



Then the expression of p e , as defined by 

Pe(x,y)-- 



1 



1 + x 2 +y 2 ' 

is used to analytically compute f e and S e with 

fe = G ePe , S e = V xPe . (4.3) 



These definitions are inserted in the numerical method described in Section [ 
to compute the numerical approximation p e ,app finally compared to the exact 
solution p e . The relative errors denoted e p , pG {2;oo}, are defined by 

\\Pe — Pe,app\\i2(j} \\Pe ~ Pe,app \\ £=o (j) 

\\Pe\\t?(I) ' V 



These quantities are displayed on Figure 1(a) as functions of the space step 
h and for different anisotropy strengths e = 10~ 1 , e = 10~ 9 and e = 0. A linear 
decrease of the errors is observed with the mesh refinement, the slope being 
equal to 2, which is consistent with the definitions (|3.2j) and (|3.3p of dh and 
dh,* as second order accurate approximations of the differential operators b • V x 
and V x -(b-). Furthermore, this property holds for all considered values of e, 
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Figure 1: Relative error \\p e — Pe,app||^-m/||Pe||.£-(r) (left) and parallel gradient 
||^/» 7r e 1 app||ft.(/ ili )/||Pe||^p(j) (right) in (i 2 and £°° norms as functions of the space 
step h for non-uniform data G £ , H e and b denned by (|4.1|) - (|4.2p . and different 
anisotropy strengths e = l, e = 10 -9 and e = 0. 



including e = 0. This demonstrates the e-invariance of the numerical scheme 
second order accuracy with respect to the space step h. 

The ability of the scheme to compute a solution component 7r e with no 
gradient in the anisotropy direction is also investigated. The numerical approx- 
imation of 7r £ , ir et app provided by (j3.4p ~ p.5p . should verify a discrete analogous 
of the property b • V x 7r e = 0. This is analyzed thanks to Figure 1(b) | where 
the evolution of ||9/ l 7r £,app||^-(j )/l|Pell^p(/) as a function of the space step is dis- 
played for p = 2,oo, e= 10 -1 , e = 10~ 9 and e = 0. Note that the quantity dhK t . ap p 
is the residual of the linear system solved to compute the solution of (|3.4|) , and 
consequently characterizes the precision of the linear system solver. For these 
test cases, a sparse direct solver being used [T], the accuracy is very close to 
the computer arithmetic precision, at least for small linear system sizes. This 
precision is observed to deteriorate moderately with the increase of the system 
size which explains the growth of the error with vanishing mesh sizes. However 
this does not affect the precision of the scheme, as demonstrated by the results 
of Figure 1(a) 



4.2 Anisotropy angle influence on the method accuracy 

In this section, we quantify the sensitivity of the numerical method with respect 
to the anisotropy direction variations. More precisely, we wish to analyze the 
accuracy of the method as a function of a, the angle measured between the 
anisotropy direction and the first direction (associated to the first coordinate). 
The anisotropy direction b is assumed to be uniform and defined as 

h= (sina, — cosa), a£[0,— ]. 
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In order to manufacture an analytic solution for the problem, we introduce 
a system of coordinates which is adapted to b. These coordinates are denoted 
(X,Y) and are deduced from (x,y) by the relations 

X = xcosa + ysma, Y = xsma — ycosa. (4-4) 

In these coordinates, the linear diffusion problem (|2.14[) writes 



-d Y (H e d Y 'P e ) + eg e Ve = eF e -dY('HeB-S e ), in Q, 
d Y V e = B-S e , ondQ, 



(4.5) 



with V e (X,Y)=p e {x,y), H e (X,Y) = H e (x,y), g e (X,Y) = G e (x,y), F e (X,Y) = 
f e (x,y), S e (X,Y) = S e (x,y) and B(X,Y) — b(x, y). It is straightforward to verify 
that the function V e given by 

V t {X, Y) = sin(X) + gX Y) dY (& Ce ) , 



is the solution of (|4.5|) provided that F e and S e satisfy 

T e = g e T £ , B-S £ = d Y V ei 

and where C e (X,Y) = l e (x,y) with l e \ gn =0. This requirement is met by the 
following definition 



le(x,y) = sin I 



2ir(x-X- 1/2 ) \ . I 27r(y-y_ 1/2 ) 



*JV e +X/2 — 2-1/2/ \ VN y + l/2 -y-1/2. 



which ensures ^(2^+1/2, yj+1/2) = for any (i,j) S I*\I*. The problem is stated 
in Cartesian coordinates thanks to the change of variables (|4.4j) yielding to 
Pe(x,y) = 7Te(x,y) + q e (x,y) with 

7r e(^ ; y) = sin(a;cosa + ysina) , q e (x,y) = — — -V x -(G e (x,y)b(x,y)l e (x,y)) , 

G € (x,y) 

the other coefficients being manufactured similarly with G e and H e given by 

In the following tests, the computation domain 0= [1,2] x [1,2] is discretized 
thanks to a uniform mesh constituted of 200 x 200 cells. The relative approxi- 
mation error as a function of the angle a is displayed on Figure [2] for different 
norms. The numerical method accuracy is observed to remain almost unaltered 
by the anisotropy direction changes. More precisely, we observe a variation of 
the relative errors in norms t 1 and I 2 lower than 4% and a variation of £°° lower 
than 7%. Furthermore, these observations are redundant for several values of 
e: in Figure [5] we have considered e= 10~ 3 and e= 10^ 8 and the obtained error 
curves are very close. Note that other experiments have been carried out for 
anisotropy strengths ranging from e = to e = 1 and with other definitions of G e 
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0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 

Angle ct 

(a) \\p e — Pe,app||#>(j)/||pe||#>(j)! e=10 3 . 




0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 

Angle ct 

(b) \\Pe ~ Pi,app\\lP{I) I llPell^P(r) i e= 10 8 - 



Figure 2: Relative error \\p e — Pe,app||^(n/]].Pe||#>(j) (p=1j2,oc) as a function 
of a: case with G e (x,y) = H e (x 7 y) = 1 + sin 2 (x) sin 2 (y) and e = 10~ 3 (left) and 
e = l(T 8 (right). 



and of H e , with comparable results. The curves being very similar to that of 
Figure [21 these plots are omitted. 

These observations confirm one of the main ideas of the present paper: the 
accuracy of the method is almost independent of the anisotropy direction rel- 
atively to the grid, i.e. the mesh over O can be constructed whatever the 
anisotropy direction and strength, without a significant loss of accuracy. 

4.3 Convergence of Gummel's loop 

The third test sequence is devoted to the convergence of the linear problems 
sequence defined in Sections 2.2 and 3.3 for solving the non-linear model (jl.ip . 
The process detailed in the preceding sections is again implemented to manu- 
facture an analytic solution for the non-linear problem. 

The computational domain remains 0= [1,2] X [1,2], the anisotropy direction 
b is a function of the space variables whose expression is given by equation (|4.2[) . 
H e and g e being defined as 

H e (x,y) = 1 + cos 2 (x) cos 2 {y), g e (p)=p 6 . (4.6) 

Note that this choice of g e (p) introduces a severe non-linearity in the prob- 
lem. Several tests have also been performed with other definitions of g e (p), for 
instance 

9e(p) = -p(l-p) (P~^j ' 

which defines an anisotropic diffusion-reaction equation similar to the steady- 
state Allen-Calm equation (see [3l [25] ) used in phase transition problems. These 
tests produce results almost identical to the results which are obtained when 
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g e (p)=p 6 so we only consider the strongly non- linear reaction term denned in 
(|4.6p within the presentation of the numerical results in the next lines. 
The solution p t is constructed thanks to a cubic spline S, precisely 

Pc( X ,y) = l+s(^^) s( V -^) , (4.7) 



L X I \ Ly 



with S(z) = for |^| ^ [0,2] and 



. i(2-|z| 3 ), if Kid < 2, , 



with (x m id,y m id) = (lit) an d = = 1/10. To analyze the convergence with 
respect to the number of Gummel's iterations, the sequence is initiated with 
Pe.o, a perturbation of the non-linear problem solution, reading 

p et o(x,y)=p e (x,y) + rjm.ax(0,l- jj,(x - x mid ) 2 - p, (y - y mi d) 2 ) , (4.9) 

where fi and rj/j, are parameters controlling the support and the magnitude of 
the perturbation. Since Gummel's method is constructed on a linearization of 
the problem its convergence cannot be guaranteed with a poor estimation of 
the solution as initial guess. It means that the parameters fi and rj cannot be 
chosen completely arbitrarily: indeed, several simulations have been performed, 
all with the same parameters except X] ranging in {0,10,20,40,60, 100, 1000} and 
fi ranging in {1, 10, 60, 100} and it has been observed that Gummel's method does 
not converge as N — >oo when r\ is larger than 10 2 . Concerning the parameter 
/i, the simulation sequence reveals that the convergence of Gummel's method is 
almost not affected by the amplitude of //. 

The successive relative errors measured between the iterates of the Gummel's 



loop and the exact solution are plotted on Figure 3(a) The computations 
are carried out on two different meshes, M 100 and M 10 oo with 100 x 100 and 
1000 x 1000 cells, with rj — 0.1 and /i = 60 and for anisotropy strengths including 
e = 0. Along with the graphical representation of the solution approximation 
error, the evolution of the corrector norm relative to that of the solution, namely 
the quantity ||5 e) i\r,applU 2 (z)/||Pe|U 2 (j) » is also plotted in Figure [3 (b)| These last 
results being almost identical for both meshes, the plot related to the finest 
mesh is omitted in this figure. 

In spite of the large perturbation amplitude, Gummel's iterative method 
converges in a small number of iterations, for both meshes and for all e-values. 
The corrector term 5 t ^N,app rapidly decreases to reach the computer precision 
threshold (10~ 15 ) after 4 iterations. In the same time, the relative error also 
decreases but the approximation is not improved by subsequent iterations, the 
error remaining constant for iteration numbers greater than 4. At this stage, 
the precision of the approximation is not limited by the linearization process of 
the Gummel's loop anymore, but by the discretization error of the linearized 
problem, explaining the plateau described by the error. To document this ana- 
lyzis further, we summarize in Tabled] the values of the relative error measured 



S. Brull, F. Deluzet, A. Mouton 



20 




Figure 3: Gummel's iteration convergence: evolution in log 10 -scales 
of the relative error E2,N=\\Pe-Pe,Nf,app\\e?(i)l\\Pe\\p{i) (left) and of 
||^e,Ar,app||£ 2 (7)/||pe|U 2 (7) (right) as functions of th iteration number N for the 
non-linear problem. The simulations are performed with the uniform meshes 
Mioo and Miqoo composed of 100 x 100 and 1000 x 1000 cells (for the corrector 
norm, the values being very similar for both meshes, only those of the coarsest 
one are displayed). 





Error 


Mioo 


M200 


M50O 


Afiooo 


io- 1 


El,N f 


3.9452 xl0~ b 


9.8116 xlO" 6 


1.5673 x 10~ b 


3.9166 x IO" 7 


E2,N f 


1.0446 xl0~ 4 


2.6188 xl0~ b 


4.1988 x 10~ b 


1.0505 xl0" b 


Eoo,N f 


6.0730 xl0~ 4 


1.5793 xl0~ 4 


2.5942 x 10~ b 


6.5451 xl0~ b 


io- 12 


El,Nf 


3.9796 xl0~ b 


9.8969 xl0~ b ' 


1.5808 x 10~ b 


3.9504 x 10" 7 


E%,Nf 


1.0496 xl0~ 4 


2.6311 xl0~ b 


4.2184 x 10~ b 


1.0554 x 10~ b 


Eoa,N f 


6.1098 xl0~ 4 


1.5885 xl0~ 4 


2.6087 x 10~ b 


6.5815 x 10~ b 





E\,N t 


3.9796 xl0~ b 


9.8969 xl0- H 


1.5808 x 10~ b 


3.9504 x 10~ 7 


E2,N f 


1.0496 xl0~ 4 


2.6311 xl0~ b 


4.2184 x 10~ b 


1.0554 x 10~ b 


Eoo,N f 


6.1098 xl0~ 4 


1.5885 xl0~ 4 


2.6087 x 10~ b 


6.5815 x 10~ b 



Table 1: Relative error E PtNf — \\p e -Pe,N f ,app\\ep(i)/\\Pe\Up(i) (p = l,2,oo) for 
the non-linear problem defined by g e (p) =p 6 . The computations are carried out 
on uniform meshes M). constituted of k x k cells (k — 100,200,500,1000) with 
several values of e and after a number of iteration of Gummel's loop Nf large 
enough for the convergence to be effective. 
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between the exact solution and the approximation obtained after Nf iterations 
of the Gummel's loop. This quantity is referred to as E p ^ s (p=l,2,oo) and 
computed for Nf large enough to ensure that the plateau above mentioned is 
reached. For the investigations carried out, this requirement is met as soon as 
Nf > 4. The approximation error E p ^ f is observed to quadratically decrease 
with the space mesh: the error norms related to the computations performed on 
a 1000 x 1000 mesh are for instance 10 2 times as small as those carried out on 
a mesh with 100 x 100 cells. This is a consequence of the second order accurate 
discretization of the spatial operator already outlined in section [4~Tl Finally, the 
results of Table Q] also demonstrate the independence of the numerical method 
precision with respect to the anisotropy intensity. 

4.4 Highlight of the scheme Asymptotic-Preserving prop- 
erty 

These last experiments are devoted to illustrate the Asymptotic-Preserving 
property of the numerical method, i.e. its ability to compute an accurate ap- 
proximation of the solution of the limit problem (|2.2|) . The solution of the 
problem is constructed as a sequence (p e ) e >o defined by 



p e =Po + ep 1 



e ' 



with 

Po(x,y) = l + sfc^j S {^j^j , (4.10) 

-1, , f n ( 2n(x-x mld ) \ ( 2n(y-y mid ) \\ 
p e (a;,y)=maxl 0,cosl — I cos I I I. (4.11) 

The functions g e , H e and b are defined as in the previous test sequence, the 
initial guess for Gummel's loop being constructed following (I4.9|) using the same 
perturbation. We now wish to evaluate the error measured between the exact 
solution of the limit problem po and the approximation computed thanks to the 
AP-schemc for vanishing e. This error, denoted E e and defined as 

E e = \\Pe,app-Po\\ i 2 (I) /\\po 



£2(J)/ 11^01^2(7) ; 



is plotted on Figure 4(a) as a function of e. The data represented on this figure 
are obtained after convergence of the Gummel's loop. Two regimes can be 
identified. The first one is related to the largest values of e for which a linear 
decrease of the error is observed. The second one is a plateau whose value 
depends on the mesh step h, this value being lower for refined meshes. Precisely 
we note a quadratic decrease of this value with the mesh size. To explain these 
features, we use the following identity 



Pe,app PO — Pe,app PQ.app + P0 
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Figure 4: Evolution in log 10 -scales of E e = \\p^ apP -f>o||.£2(j) / lboll £ 2 (/) (left) and 
of E e , app = \\pe,a PP -Po, appWpin /lbo||^(/) (right) as functions of e computed 
on uniform meshes constituted of 200 x 200 cells (M 2 oo) and 1000 x 1000 cells 
(Af 10 oo) are considered. The simulations are performed with b, H e , g e defined 
by d^B]) and 



This yields E e < E^ app + e where e = ||po,a PP -Po||p (/) / |bo||^ (/) represents the 
approximation error of po, Po : app being the numerical approximation of po P r0 ~ 
vided by the AP-scheme with e = 0, and E^ app = \\p e ,app -Po,app\\p{i) I lboll^(/)- 
The error E e linearly decreases with e as long as the approximation error eg is 
negligible compared to E e ^ app (see Figure [4(b)] ) . Below a given e- value, varying 
with the mesh size, the total error can be assimilated to eo and the decrease 
of e is ineffective. The discrete operators being second order accurate eo is 
quadratically decreasing with the mesh step h. 

As a consequence, we can conclude that p e ,app converges to po when e con- 
verges to alongwith h. This is exactly the Asymptotic-Preserving property of 
the scheme we intended to validate. 



5 Conclusions and perspectives 

In this paper we have presented an Asymptotic-Preserving numerical method for 
singular perturbation of non-linear anisotropic reaction-diffusion problems. The 
Asymptotic-Preserving property of the scheme is ensured thanks to a solution 
decomposition explained in full details in the most simple framework of a linear 
problem. This method is then generalized to non-linear problems thanks to 
Gummel's linearization method. 

In a second part, several two-dimensional numerical investigations of the 
AP-scheme are performed. These tests reveal a very weak dependence of the 
scheme accuracy with respect to the anisotropy direction, demonstrating the 
relevance of the use of non-adapted coordinates. The Asymptotic-Preserving 
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property of the scheme is also validated for vanishing e on linear as well as 
non-linear problems. The solution of the limit problem is accurately captured 
with no restrictions on the anisotropy strength. Furthermore, the computational 
efficiency of the method, in terms of memory as well as CPU usage, does not 
depend on this anisotropy strength. 

Several applications of the present work can be investigated: at present 
time, the method has been used for the resolution of linear anisotropic diffusion 
problems for a two-fluid Euler-Lorentz model (see [6]) and the non-linear version 
of the method will be coupled to an Asymptotic-Preserving scheme for a one- 
fluid full Euler-Lorentz model (see [18]). 
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